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We study the distribution of the time-integrated current in an exactly-solvable toy model of heat 
conduction, both analytically and numerically. The simplicity of the model allows us to derive the 
full current large deviation function and the system statistics during a large deviation event. In 
this way we unveil a relation between system statistics at the end of a large deviation event and 
$H , for intermediate times. Midtime statistics is independent of the sign of the current, a reflection of 

the time-reversal symmetry of microscopic dynamics, while endtime statistics do depend on the 
current sign, and also on its microscopic definition. We compare our exact results with simulations 
based on the direct evaluation of large deviation functions, analyzing the finite-size corrections of 
1^-^ I this simulation method and deriving detailed bounds for its applicability. We also show how the 

Gallavotti-Cohen fluctuation theorem can be used to determine the range of validity of simulation 
results. 

o 

S" 

'. I. INTRODUCTION 
a ■ 

C/3 , Many nonequilibrium systems typically exhibit currents of different observables (e.g., mass, energy, spin, etc.) 
■ which characterize their macroscopic behavior. These currents fluctuate, and the (large deviation) functional which 
determines the probability of these fluctuations is a natural candidate to generalize the concept of free energy to 
^ nonequilibrium systems. In this way, understanding how microscopic dynamics determine the long-time averages of 
these currents and their fluctuations is one of the main objectives of nonequilibrium statistical physics ([it 12; [SyJ; 
H; S 0; S S EO)- An important step in this direction has been the development of fluctuation theorems ((a: lol: liol) . 
Q • which relate the probability of forward and backward currents reflecting the time-reversal symmetry of microscopic 
Q \ dynamics. These are relations of the form 

I where P{q,t) is the probability of observing a time-integrated current Qt = after a long time t, and E is some 
. constant such that Eq corresponds to the rate of entropy production (0; S H [13; [HI) ■ Despite this progress, we still 
lO ' lack a general approach based on few simple principles to understand current statistics in nonequilibrium systems. 

This has triggered an intense research effort in recent years. For instance, Bertini, De Sole, Gabrielli, Jona-Lasinio and 
Landim (0; y) have recently introduced a hydrodynamic fluctuation theory (HFT) to study both current and profile 
00 i large dynamic fluctuations in nonequilibrium steady states, providing a variational principle which describes the most 
' probable (possibly time-dependent) profile responsible of a given current fluctuation. Simultaneously, Bodineau and 
^ [ Derrida (jl|) have conjectured an additivity principle for current fluctuations which enables one to explicitely calculate 
I the current distribution for ID diffusive systems in contact with two boundary baths. This principle, which is 
p\ ' equivalent within HFT to the h ypo thesis that the optimal profile is time-independent (Q), has been recently confirmed 
in a model of heat conduction (l2). 

The predictions derived from the additivity principle and the hydrodynamic fluctuation theory have been tested 
in few models for which the exact current distribution can be calculated. For more general cases in which no exact 
solutions are available, new algorithms that allow the direct evaluation of large deviation functions in computer 
simulations have been recently proposed (flst These methods are based on a modification of the microscopic 

dynamics, so that the rare events responsible of the current large deviation are no longer rare. Though successful, the 
new algorithms seem to suffer from finite-size effects when measuring the probability of extreme current deviations 
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(fl3 ). Here again the existence of simplified models for which exact solutions can be obtained is essential to elucidate 
the origin and importance of these finite-size corrections in simulations. 

The aim of this paper is to present one of these simplified models, i.e. a toy model of heat transport between two 
thermal baths at different temperatures. The simplicity of this model allows us to gain a better understanding of: (i) 
system statistics during a large deviation event and how it reflects the time reversibility of microscopic dynamics, and 
(ii) finite-size corrections to the direct evaluation of large deviation functions. The model and the calculation of its 
current large deviation function are presented is Section II. The probabilities of having certain configuration at the 
end of a large deviation event and for intermediate times are calculated in Section III, while Section IV is devoted to 
analyze their relation in terms of the reversibility of microscopic dynamics. Simulation results are described in Section 
V, together with an analysis of finite-size corrections affecting the direct evaluation of large deviation functions. This 
analysis provides a deep understanding of the range of validity of this simulation method, which otherwise can be 
determined using the Gallavotti-Cohen fluctuation theorem. Finally, Section VI contains our conclusions. 



II. MODEL AND EXACT SOLUTION 

Our toy model consists in a single lattice site characterized by an energy e G IR+. This site is coupled to two 
thermal baths at different temperatures, (left) and Tr (right), and we assume > Tji without loss of generality. 
Dynamics is stochastic and proceeds as follows: (i) With equal probability, we randomly choose the system to interact 
with one of the heat baths, Tl or Tr; (ii) a new energy e' is drawn from the distribution f3i,ex.p{—P^e') corresponding 
to the (inverse) temperature of the selected heat bath, /3y = T~^, with i/ = L,R. The presence of a temperature 
gradient then forces the system out of equilibrium. The transition rate U^l'^ from configuration e to e' interacting 
with bath v = L, R can be written as 

and Ue'e = C^e'^^ + ^e'e^ ■ ^hc transition rate is normalized, J2e' ^e'e — 1, thus guaranteeing the conservation of prob- 
ability during the stochastic evolution. This model can be regarded as the single-site version of the one-dimensional 
Kipnis-Marchioro-Prcsutti (KMP) model of heat conduction (12; 15). Having a single site then implies that no energy 
diffusion takes place, but only energy exchange with the thermal reservoirs. Associated to the interaction with the 
baths there is an energy current through the system, g^^g. We may define this current at the microscopic level in 
different ways. For the time being, let us define q^^fl in a symmetric way 




In this way, positive currents correspond to the transport of energy from the left to the right reservoir. In the steady 
state, the probability density for having an energy e can be trivially calculated 

P(e)-i(/3ie-^-^-H/3fle-'3-^) , 

and hence local equilibrium does not hold for this toy model. Notice that P(e), despite being non-Gibbsian, can be 
related to the transition rate Ue'e via detailed balance. The mean energy is (e) = ^C^l + Tji), while the average 
current is (q) = j{T]^ — 7^). Therefore Fourier's law trivially holds in this system for arbitrarily large temperature 
gradients, VT = ^(T^ — T^), with a thermal conductivity k{T) = i. Moreover, in equilibrium {Tl = = T), 
current fluctuations have a variance (t{T) = T^. These values of k{T) and cr(T) agree with the thermal conductivity 
and current variance of the spatially-extended KMP model ([ish . reinforcing the view that this model is just the single 
site version of the KMP process. 



A. Current Large Deviation Function 



Let now P(e, q, t; cq) be the probability of having the system in configuration e with a total time-integrated current 
qt after a long time t, starting from a configuration cq. This probability depends also on the temperatures of the 
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boundary thermal baths, but we drop this dependence here for notation convenience. P(e, q, t; Bq) obeys the following 
master equation^ 

P(e, q, t; eo) ^Y.T. ^eM^' , 9* - , ^ " 1; eo) • 

e' I' 

Iterating in time, we can write the above probability as 

P(e, q, t- eo) - E • • ■ f^i'o - i^'feL + • ■ • + • 

et_i...ei vt-.-vi 

This is nothing but the weighted sum over all possible phase-space paths {e . . . eg} starting at eg and ending at e, 
with duration t, such that the total time-integrated current is qt. One may then write the probability of obser ving a 
total current qt as P(g, t; ep) — 9' ^5 ^o)- For long times P(g, t; Bq) obeys a large deviation principle ([l6l : flTl ) 

P(<Z,t;eo)~e+*^(^), (3) 

where is the current large-deviation function (LDF), which does not depend on the initial state bq. The function 
J'{q) is typically everywhere negative except for q = (q), for which it is zero, meaning that current fluctuations away 
from the average are exponentially unlikely in time. In most cases it is convenient to work with the moment-generating 
functions of the above distributions 

Uie,X,t;eo)^Y.'^''^P(^'q,t;e,)= Y ^^ll ■ ■ ■ ^^^1 ' (4) 

q et-i...ei Vt...vi 



and n(A,t;eo) = n(e. A, t; eo). For long t, one can show that n(A,t;eo) ^ e*^^'^\ where ^(A) = maxg[J^((7) + \q] 
is the Legendre transform of the current LDF. We can now define a modified dynamics, U^^l = e'^'^<:'= u'^^l, so 



e...ei i^f-i^l 

This dynamics is however not normalized. For our particular model 

r ^e-('^-^)^'e-^^ ,. = L 
u'jI-{ ' (6) 

The simplicity of this single-site model allows the direct calculation of the sum ([5]). Let us start by defining the 
following recurrence 

et-i vt-l 

with xj'^-'(e; eo) = uiel- In this way, eq. ([5]) can be written as 



n(A,t;eo)=5]^x(^'(e;eo). (7) 

e y 

A simple induction argument shows that 



All throughout the paper we are using the symbol J^^, in a general context, meaning "sum over configurations" , independently of whether 
configuration space is discrete or continuous. In fact, for the model studied here configuration space is continuous, and sums become 
integrals over the system energy. 
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where A[^'^\eo) are some coefficients, and —(3r < X < Pl for the integrals to converge. Plugging this into into eq. 
([7]) we arrive at 



n(A,t;eo) 



{L)i 



The coefficients A\^^ (eo) obey the following recurrence (in matrix form) 



At+i 




1 

V /3fl + A 



1 



(8) 



A\ 



with ^^■'^■'(eo) — f3j^e2'^°/2 and A[^\eQ) = P^e 2'^"/2. The recurrence matrix S has eigenvalues a± = [1 ± 
(/^H+Alffe-A) ]/^' ^ ^^^^ associated eigenvectors ~ [1,±(/)(A)], with 



V /3ii(/3fl + A) 



Using these eigenvectors as a basis, we write ^i(eo) = a+(eo)V^+ + a-(eo)V'-, so = S** -^Ai 

a_Q!^^'0_. For long times the largest eigenvalue dominates, so a\.^^ ~ a+a^"'^ and A^'"^ ^ a-\-(j){\)a^^ in this 
limit. Using these long time asymptotics in eq. ([5]), and recalling that n(A,t; eo) ~ exp[t/i(A)] for long times, we find 



/i(A) = In . 




{(iR + \){f3L-\) 



(9) 



Fig. [T]shows /i(A) versus A//?/? for different values of the ratio Pl/ Pr.- For all pairs we have /x(0) = due 

to the normalization of P(q, t; eg), see the definition of n(A, t\ eo) in eq. (H]) and discussion below. One can also check 
that the first derivative of /z(A) evaluated at A = yields the average current (flTt). i.e. /x'(0) = {q) = \{Tl — Tr). 
In addition, as a result of the time reversibility of microscopic dynamics, the Gallavotti- Cohen fluctuation relation 
([sl: [9I: [loi) holds for this system, /i(A) = iJ,{—X — E) with iJ = Pr—_0l- This relation for /i(A) is equivalent to the usual 
fluctuation relation J^(-g) = J^{q) - Eq for the current LDF |3; SB [l3). 
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FIG. 2 Ratio of average energy at the end of a large deviation event and for intermediate times, for different values of 
A e and [3l G [0, 1], with fixed Pr = 1. 



III. STATISTICS FOR A LARGE DEVIATION EVENT 



We now study the system energy distribution in a large deviation event of (long) duration t and time-integrated 
current qt. One may measure this distribution at a time r such that: (i) 1 <gC r <C i.e. for intermediate times, or 
(ii) for T ~ t, i.e. at the end of the large deviation event. We first focus on endtime statistics. Consider then the 
probability Pend(e|g; t, eo) = P(e, g, eo)/P((j', t; eg) of having a configuration e at the end of a large deviation event 
associated to a current q. One can easily show in the long time limit that Pcnd(e|A; t, eo) = n(e, X,t; eo)/n(A, t; Cq) = 
Pend[e|(7*(A); i, eo], where q*{X) is the current conjugate to parameter A, such that J-'{q*) + A = 0, with n(A,<;eo) 



defined in eq. ([5|) and n(e. A, t; eo) = '{e; eo)+Xj: '{e; eg). For long times Pcnd(e|A; t, eg) converges to the following 
distribution, independent of both t and the initial state eq. 



Pond(e|A) = Rx 



0(A)e 



(10) 



where i? is a normalization constant, and we have used that a[^^/a[^^ for t > 1. 

To compute the energy distribution for intermediate times, 1 t <^ t, let P{e, q,T,t; cq) be the probability that 
the system was in configuration e at time r when at time t the total current is qt, starting from configuration eg. This 
probability can be written as 



P{e,q,T,t; eo) 



E E u^t. ■ ■ ■ utvJ uizl ■ ■ ■ u^e:il s[qt - {qi:ii^ + ...+ qtl)] 



et...e-r+ie-r-i...ei Ut...vi 



where we do not sum over e. Defining the moment-generating function of the above distribution, n(e. A, r, t; eo) = 
exp(iAg)P(e, g, T, eo), we can again check that the probability weight of configuration e at intermediate 
time r during a large deviation event of current q, Pi„id(e|q; r, eo) = -P(e, g, r, t; eo)/P(q, i; gq), is also given by 
II(e, A, r, i; eo)/n(A, eo) for long times such that 1 ^ t ^ t, with q — q*{X)- Now, 



n(e,A,r, i;eo) 



eo 



(11) 



so we can write Pmid(e| A; t, t, eo) = Il{\,t — r; e) x n(e. A, r; eo)/n(A, t; eo). Therefore, using eq. jS]) in the limit of 
long times t and r, such that 1 ^ r ^ t, we arrive at 



Pmid(elA) = i?' X 



e^'^^" + 0(A)e-('^^-^)'= + (/.(-A - £;)e-('3«+^)" + ^e'^"- 



(12) 
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which again does not depend on Cq and on times t and r (as far as r ^ t) . Here R' is another normalization constant 

e~i'^). Using eqs. PH)) and p^ . we obtain the system average 



and we have used the fact that a+(e) — (^ea' 



40(A) 



energy at the end of the large deviation event and for intermediate times 

A 



(e)ond(A) = 



(e}mid(A) = 



A 



A. 



Tr 



4a+(A) 



I3r 



A 



) + 0(A)(/3j^ + ^) 



0(A) 



Pi 



Tl 



(13) 



(14) 



(/3l - A)2 ' ' {13r + A)2 0(A) 

Fig. [2] shows the ratio [(e)cnd/ (e)mid](A) as a function of A and /3l, for /3fl = 1. This ratio is typically smaller than 1, 
except for A « (i.e. small fluctuations around the average current). Therefore the average system energy at the end 
of a large deviation event systematically underestimates the average energy for intermediate times. This results from 
the slower exponential decay of Pmid(e|A) for large energies as compared to Pcnd(e|A). 

We could do again all calculations using a different microscopic definition of the current q^if,, see eq. ([2]). As far as 
the internal energy of the system remains bounded, the difference between the time-integrated currents obtained with 
different definitions of qe'e will be also bounded, thus guaranteeing the same long-time average current (0)- Therefore 
any observable obeying the Gallavotti-Cohen symmetry should be independent of the microscopic definition of the 
current. For instance, defining the current as the energy exchanged with the left heat bath 







Left bath 



Right bath , 



(15) 



we obtain the same results for the large deviation function, eq. ([9|), as expected. Moreover, results for the midtime 
energy statistics do not change for the new current definition. However, results at the end of the large deviation event 
do change with respect to the symmetric current definition, eq. ([2|). In particular, we obtain for the left current (jlSp 



end 



(e|A) = i?x 



(e}cnd(A) = 



A. 
2' 



A 



) + 0(A)/3fl 



which should be compared with eqs. (HU]) and p^ . respectively. This sensitivity of endtime statistics to the micro- 
scopic definition of the current is a reflection of the violation of the Gallavotti-Cohen symmetry for this distribution 
(see below). 



IV. TIME REVERSIBILITY AND RELATION BETWEEN Pmid(e|A) AND Pe„d(e|A) 

Direct inspection of eq. reveals that Pmid(e|A) = Pmid(e| — A — _E), with E — Pr — Pli or equivalently 

Pmid(e|(7) — Pmid(e| ~ ?), SO systcm statistics at intermediate times during a large deviation event does not depend on 
the sign of the current. This implies in particular that (e")niid(9) = (e")mid(— </) Vn. This counterintuitive symmetry of 
Pmid(e|9) is a reflection of the time reversibility of microscopic dynamics, and it's a general result for systems obeying 
the local detailed balance condition (fH). For our particular system this condition reads C/ee'Peq(e') = Ue' ePeci_{e)e~^'^'' ' , 
where Pcc^e.) = ^-"^^-^ exp[— ^-"+^^ e] is an effective equilibrium weight for configuration e. Local detailed balance then 
implies a general symmetry between the forward modified dynamics for a current fluctuation and the time-reversed 
modified dynamics for the negative current fluctuation, i.e. L'ee'(A) = P(Jq^(e')J7e'e(— A — E)pcq{e), or in matrix form 

[/T(A)=P-ql[/(-A-i?)Peq, (16) 

where U"^ is the transpose of U and Peq is a diagonal matrix with components Peq(e) pl[ ). Introducing 
Dirac's notation^ we can write the spectral decomposition of the modified dynamics (flit [l3t [ist [T9I ). U{X) = 



^ Dirac's bra and ket notation is useful in the context of the quantum Hamiltonian formalism for the master equation Il3t llSt Il9l) . 

The idea is to assign to each system configuration e a vector [e) in phase space, which together with its transposed vector (e|, form 
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e^3(-'*)|Aj^(A))(A^(A)|, where we assume that a complete biorthogonal basis of right and left eigenvectors for 
matrix U{X) exists, U\Af{X)) = e^^^^^\Af{X)) and {Af{X)\U = c^^^^'> {Af{X)\. Denoting e^(-^) the largest eigenvalue 
of {/(A), with associated right and left eigenvectors |A^(A)) and (A^(A)|, respectively, we find for long times 

n(e, A,t;eo) = {e\U*\eo) ^ e*^W(A^(A)|eo>(e|A«(A)) . (17) 

In this way we have fJ.{X) = A(A), i.e. the Legendre transform of the current LDF is given by the natural logarithm of 
the largest eigenvalue of U{X). Since eq. p6|) implies that all eigenvalues of U{X) and [/(—A — E) are equal, and in 
particular the largest, we find that /x(A) = A — E) and this proves the Gallavotti-Cohen fluctuation relation. Now, 
eq. ([T7|) also implies that Pend(e|A) cx (e|A^(A)) for long times, so the right eigenvector |A^(A)) associated to the 
largest eigenvalue of matrix U{X) gives the probability of having a configuration e at the end of the large deviation 
event. In a similar way, one can easily show from eq. (jlip that, in the long time limit, 

P,„id(e|A)oc(A^(A)|e)(e|A^(A)}. 

In general, the eigenvectors |A''^(A)) and (A^(A)| are different because U is not symmetric. In order to compute the 
left eigenvector, notice that |A^(A)) is the right eigenvector of the transpose matrix J7"^(A) with eigenvalue e^^'*''. 
Therefore, using the symmetry relation eq. p6|) . one can show that if |A^(— A ~ E)) is the right eigenvector of 
f7(-A - E) with largest eigenvalue, such that |A^(-A - E)) = Ee(e|A^(-A - E))\e), then 

\A^{X))=J2p-^\e){e\A^i-X-E))\e) 

e 

is the right eigenvector of U^{X) associated to the same eigenvalue. In this way we find 

P.„M(e|A) (X (A^(A)|e)(e|A«(A)) ^ p;^\e){e\A^{-X - E)){e\A''iX)) , 

or equivalently 

P.M(e|A) = K X Pend(e|A)P.„,(e|-A^i;) 

Pcq(e) 

with K some normalization constant. Hence important configurations at intermediate times are those with an sig- 
nificant probabilistic weight at the end of both the large deviation event and its time-reversed process. The above 
equation relates midtime and endtime statistics, and proves that in general Pmid(e|A) = Pmid(e| — X — E). 

Let us end this section by pointing out that the computational method of Ref. (fisl). which allows the direct 
evaluation of large deviation functions in simulations (see next section) , yields also the endtime statistics associated 
to a large deviation event, but this method cannot access midtime statistics (see however Ref. (|20l )). In this way, eq. 
(IT5| can be used in simulations to obtain midtime statistics from endtime histograms, as done below. 

V. FINITE SIZE EFFECTS DURING THE DIRECT EVALUATION OF LARGE DEVIATION FUNCTIONS 

In general, measuring in simulations the current large deviation function (or any other type of LDF) is a very 
difficult task. This is because, by definition, LDF's involve exponentially unlikely events, see eq. To overcome 
this problem, Giardina, Kurchan and Peliti recently proposed an efficient computational scheme to measure LDF's 
in many particle systems (fisl ). This method uses the modified dynamics Ue'e defined in eqs. ©-([G]), for which the 
rare events responsible of the large deviation of the observable are no longer rare (fisl : IT3 ). Despite its success, this 
simulation method fails to provide accurate measurements of the LDF for extreme current fluctuations (fl^. In this 
section we show, using our toy model as an illustrative example, that these deviations are generally present, and can 
be traced back to finite-size effects associated to the direct evaluation of large deviation functions. Furthermore, the 



an orthogonal basis of a complex space and its dual jlSt [l9h . For instance, if the number of configurations available to our system 
is finite (which is not the case in this paper), one could write |e)^ = (e| = {. . . . . . 0, 1, . . . . . .), i.e. all components equal to zero 
except for the component corresponding to configuration e, which is 1. In this notation, U^/g, = {e'\U\e), and a probability distribution 
can be written as a probability vector \P{t)) = '^g, P{e,t)\e), where P{e,t) = (e|P(t)} with the scalar product (e'|e) = S{e' — e). If 
(s| = (1 . . . 1), normalization then implies {s\P(t)) = 1. 
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Mt 



FIG. 3 Sketch of the evolution and cloning of the copies during the direct evaluation of the large deviation function. 



simplicity of our toy transport model allows for an accurate estimation of the range of validity of the method of Ref. 
(fish, and the scaling of this range with the size parameter. We also show that the Gallavotti-Cohen symmetry can 
be used to bound numerically this range of validity. 

The direct evaluation of large deviation functions is based on the sum 

et...ei 

where TJe'e = U^f^ , see eq. (HI), and tjj:^^ — uj,^^ e:xjp[Xql!fl] is the modified dynamics. In principle this sum over 

paths could be easily computed by a Monte Carlo sampling with transition rates Ue'e- However, the dynamics U is 
not normalized, so we introduce the exit rates Ye — X^e' X^i/ ^e'^e define the normalized dynamics U^t^ = Y^^tJe'e- 
With this definition 

n(A,i;eo)= Ye,-iUU^,---YeoK,eo- (19) 

et...ei 

This sum over paths can be realized by considering an ensemble of M 3> 1 copies (or clones) of the system, evolving 
sequentially according to the following Monte Carlo scheme (fisl : 

I Each copy evolves independently according to modified normalized dynamics U' . 

II Each copy m G (in configuration et[m] at time t) is cloned with rate Yej[m]- This means that, for each 

copy m, we generate a number ifet[m] = L^et[?n]J + 1 of identical clones with probability Yetfrix] ~ L^et[r?x]Ji or 
^etlm] — L^et[m]J Otherwise (here [xj represents the integer part of x). Note that if -ftTetfrn] — the copy may 
be killed and leave no offspring. This procedure gives rise to a total of M/ = J2m=i ^et[m] copies after cloning 
all of the original M copies. 

Ill Once all copies evolve and clone, the total number of copies Mj. is sent back to M by an uniform cloning 
probability Xt = M/M[. 

This process is sketched in Fig. [3l It then can be shown that, for long times, we recover /Lt(A) via 

Ai(A) = -iln(Xf •••Xo) fort>l (20) 

To derive this expression, first consider the cloning dynamics above, but without keeping the total number of clones 
constant, i.e. forgetting about step III. In this case, for a given history {et, et-i • . . ei, cq}, the number M{et • . . cq; t) 
of copies in configuration at time t obeys M(et . . . Cq; i) = Ul^^^__Ye^__^M {et-i ... eg; t — 1), so that 

M{et ...eo;t)= r,,„,[/;,^_^ . . . r,„C/;,„AA(eo, 0) . 
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FIG. 4 Left panel: Simulation results for ^i(A) and the exact curve. The inset shows the difference between theory and numerical 
results. Right panel: We test the Gallavotti-Cohen relation by plotting together /i(A) and A — E), with E — I3r ~ Pl- The 
inset shows their difference. 




FIG. 5 Left: Energy histograms measured at the end of the large deviation event for different values of A. Points are simulation 
results, and lines are exact results, eq. (|IOp . Right: Energy histograms at intermediate times, obtained from endtime statistics 
via eq. (|I8p . Lines are exact results from eq. (|12|l . 



Summing over all histories of duration t starting at eo, see eq. (|19|) . we find that the average total number of clones 
at long times shows exponential behavior, {Af{t)) = J^et ei -^(e* . . . eo; i) ^ M{eo, 0) exp[tfi{X)]. Now, going back to 
step III above, when the fixed number of copies M is large enough, we have Xt — {J^{t — l))/{J\f{t)) for the global 
cloning factors, so Xt ■ ■ ■ Xi = JV{eo,0)/ {N'{t)) and we recover expression pO| for /i(A). In addition, the fraction 
of clones in configuration e among the M existing copies corresponds, for large M, to the ratio {J\f{e,t)) / {J\f{t)), 
where {Af{e, t)) = X]et_i a -^(^j ^t-i ■ ■ • eo; t), and this ratio is nothing but Pend(e|A), see discussion before eq. ((TU]). 
Therefore, the fraction of clones in a given configuration yields the system endtime statistics. Finally, using the 
relation between Pmid(e|A) and Pcnd(e|A), we may also measure midtime statistics. 

We used the computational scheme described above to obtain numerical estimates of /i(A) and both Pmid(e|A) and 
Pond(e|A) for our toy model. In particular, we performed simulations for T/, = 2 and — 1, using the symmetric 
definition of the microscopic current, see eq. ([2]). Left panel in Fig. [4] shows a comparison between simulation 
results for /i(A) and eq. Numerical measurements nicely reproduce the exact result, except for very large current 
deviations, see inset to Fig. |4] (left). Note that these slight differences for extreme events seem to occur earlier for 
currents against the gradient, i.e. for A < 0. On the other hand, Fig. [5] shows the energy histogram measured at 



1,8 



S1,6 

A 
CD 
V 



1,4 



1,2 



1 1 1 1 1 1 1 1 1 
(T +T)/2 


' 1 ' 1 





1 . 1 



1 -0,8 -0,6 -0,4 -0,2 



0,2 0,4 



— Exact <e>^ 

— Exact <e> 




"1 ' r 



end 

mid 



O <e>mid[-^-(PR-PL)] 



-1 -0,8 -0,6 



-0,4 -0,2 



J 1 L 



0,2 0,4 



FIG. 6 Left: Average energy measured at the end of the large deviation event as a function of A. Lines are exact predictions, 
see eq. (|13[l . Right: Average energy at intermediate times, obtained from endtime statistics via eq. (|18|l . compared with the 
exact result, eq. (fT4)) . (e)inid(— A — E), with E = (3r — f3L, is also plotted to check the symmetry of the average energy at 
intermediate times. For completeness, (e)cnd(A) is also shown. 



the end of the large deviation event (left) and for intermediate times (right), for different values of A, while Fig. [6] 
shows the average endtime (left) and midtime (right) energies. In all cases the agreement with exact results is again 
excellent, and deviations from the theoretical predictions are only apparent for extreme values of A. Our aim in what 
follows is to show that these deviations are due to finite size effects in the simulation method described above (isl). 

The direct evaluation of the large deviation function is exact in the limit A/ — » oo. For a large but finite number of 
copies M, this method can fail whenever the largest exit rate yet[m] among the set of M copies becomes of the order 
of M itself. This condition implies that, after the cloning procedure (see Fig. [3]), configuration et[m] will overpopulate 
all the other copies, hence introducing a bias in the Monte Carlo sampling. We can estimate at which point this source 
of error becomes dominant by using extreme value statistics (j2ll ). In this way, let y^™"^ = max(Ye4[i], • ■ • ,Yet[M]) be 
the maximum exit rate among the set of M copies at a given time. The probability that is smaller than a given 

threshold y can be written as P^'^^(y, A) — [1 — P>(7/, A)]^^, where P>(y, A) is the probability that the exit rate of 
a copy is larger than y, and here we assume statistical independence among the different copies. For large M, the 
tail of P>(y, A) dominates [i.e. small values of P>(t/, A)], and P™"'^(y, A) ~ exp[— MP>(j/, A)]. Therefore, the value 
y*{X,p) of the maximum which will not be exceeded with probability p satisfies the following equation 

P>r(A,p),A] = i-ln(i). 

M p 

Now, the maximum exit rate allowed by the algorithm is M itself. Setting y*{\c,p) = M in the previous identity, 
we get an equation for Xc{p), the critical value of A beyond which the maximum which will not be exceeded with 
probability p is larger than AI, 

P>[Af,A,(p)] = i-ln(i). (21) 
M p 

This condition signals (with confidence level p) the onset of the systematic bias due to the finite number of clones 
in simulations. To further proceed, we need the probability P>(j/, A) that the exit rate of a copy is larger than y. 
For that, we first recall that the fraction of copies in a state e is just {Af{e,t)) / {Af{t)) = Pcnd(e|A), see Section III. 
Therefore the probability density for having an exit rate — U^'e during the numerical evaluation of /i(A) is 
P(Fe,A) = Ee" Pendie"\X) S[Y,-Ye>>], and 



P>(y,A)= / dyeP(i;,A) = ^Pend(e"|A)(?(i;" -y) 



where 0{x) is the Heaviside step function. Therefore the knowledge of endtime statistics during a large deviation 
event allows for an estimation of the range of validity of the simulation method of Ref. (flsl V 
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FIG. 7 Left; Simulation results for ^t(A) obtained for different number of clones M, and exact result. The inset shows the 
difference between theory and numerical results for different M. Right: Test of the Gallavotti-Cohen symmetry for different 
values of M. 



We can make these calculations explicitely in our toy transport model, for which we know Pcnd(e|A), see eq. (|10l) . 
From the definition of Ue'e for our model, eq. ([6]), we find 



2(3l-\ 



e"2<= + 



Solving for e in this equation 



2(iR + \ 



-e2- 



e±(i;,A) = -ln 



fiR + A/2 

Pr 



PrPl 



(/3fl, + A/2)(/?i-A/2) 



so in principle only two possible energy configurations, e±, are compatible with a given exit rate Y^. For large 
values of Y^., which correspond to the dominating tail of ^(Y^, A), there is only one solution with physical meaning 
(i.e. positive energy): e^iY^.X) for A > and e_(ye,A) for A < 0. Therefore, by conservation of probability, 
P(Ye, A) dl^e = Pond (e I A) de, and in the limit of large Y^ we find the following asymptotic scaling 



Pr 



p(re,A>o)-r; 



< A < 



p(i;,A<o)' 



Ye 



- 2|A 



-/3«<A<-^, 



(22) 



These power laws imply also algebraic behavior for the cumulative distribution, P>{y, X) ^ y '"'^'^^ ^1 if a{X) is the 

ln[ln(p-i)] 
In(Af) ' 



exponent for P(Fe, A). Therefore from eq. (|2T|) we find a{Xc) = 2— , arriving at'^ 



AJ(P) 



21n[ln(p-i)] 
In(Af) 



Pr 



1 21n[ln(p-i)] 
In(A.f) 



-Pr- 



(23) 



In this way, for a simulation with M copies and A G A^ (p), A^(p)], the maximum exit rate among all the copies 
will be smaller than M with a probability p, thus guaranteeing that (with confidence level p) no finite size effects 



^ This calculation neglects the A-dependent amplitudes of the power-laws P{Ye, A) in eq. I l22t . These constants give subleading corrections 
to the estimation of A^(p) derived from the A-dependence of the power-law exponents. On the other hand, notice that in order to obtain 
Xc (p) we use a{X) = because the alternative exponent «(A) = ^^^u i '^^ yields values of AJ(p) outside its domain, — < A < 0. 
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will bias the result. Therefore the interval [— A~(p), A+(p)] defines the range of validity of the method of Ref. (fT^ . 
For a high confidence level p — 0.99 and M — 10"^ clones (as used in our simulations), we find X^ip) ~ 0.38 and 
K-{P) ~ —0.75, in excellent agreement with the observed deviation from the exact result, see inset to left panel in 
Fig. m This calculation can be repeated for different values of Af, finding equally good agreement as far as M is large 
enough, see Fig. [7] (left). 

An important consequence of eqs. (l23|) is that the range of validity of the algorithm, [— A^(p), A+(p)], increases 
logarithmically with the number of clones used in the simulation. Therefore an appreciable increase of this range of 
validity in A-space demands an exponential increase in the number of clones, which is unfeasible in most situations of 
interest. We conjecture that this logarithmic increase of the range [— Aj(p), A+(p)] is not particular of our toy model, 
but a general feature for models with bounded values of A. However, this shortcoming does not forbid in practice 
to push the algorithm far enough to observe very large current deviations. Indeed, the current q*{\) conjugate to 
parameter A increases very rapidly as A approaches its bounds (-/3_r and (3l in this model), so a small (logarithmic) 
increase of the validity range [— Aj(p), A+(p)] for /i(A) is translated into a large increase of the range of current values 
for which we can measure with confidence the current LDF, J-{q). 

For most many particle non-trivial systems of interest we usually don't know the analytical form of the large 
deviation function to be measured in simulations. It is hence necessary to develop a reliable method to determine the 
range of validity of our numerical measurements. It is at this point where the symmetries of the large deviation function 
play a prominent role. In particular, most systems of interest obey the Gallavotti-Cohen fluctuation theorem which 
results from the time reversibility of microscopic dynamics. This relation can be stated as /i(A) — /i(— A — E) for the 
Legendre transform of the LDF, with E some constant, see also eq. (H]). For our toy model, E — Pr — /3l- Violations of 
this fluctuation symmetry are expected whenever finite-size corrections induce a bias in the measurement of /^(A), and 
therefore the Gallavotti-Cohen symmetry of the current LDF can be used to numerically bound the range of validity 
of simulation results. This is confirmed in Figs. 2] and [7] (right panels), to be compared with their corresponding left 
panels. 

VI. CONCLUSIONS 

In this paper we studied a toy model of heat transport between two thermal baths at different temperatures. The 
model is simple enough so the exact current large deviation function can be analytically calculated, as well as the 
full system statistics associated to a large deviation event. In this way we find a relation between system statistics 
at the end of the large deviation event and for intermediate times, which shows that important configurations at 
intermediate times are those with a significant probabilistic weight at the end of both the large deviation event and 
its time-reversed process. Thus midtime statistics turns out to be independent of the sign of the current, a reflection 
of the time-reversal symmetry of microscopic dynamics, while endtime statistics does depend on the current sign, 
and also on the microscopic definition of the current. The relation between midtime and endtime statistics is a 
general property for systems obeying the local detailed balance condition, which guarantees the time-reversibility of 
the dynamics. We also compared our exact results with simulation data to analyze the range of validity of a recently 
proposed computational scheme to directly evaluate large deviation functions. This comparison offers insights into 
the finite-size corrections associated to this simulation method. In particular, we were able to calculate the range 
of validity of the numerical results for our simplified model, finding that this range grows logarithmically with the 
number of clones involved in the evaluation of the current large deviation function. We conjecture that such slow 
growth, which impedes in practice obtaining reliable results for extreme current fluctuations, is a general feature of 
the simulation method. Finally, for more general systems for which we do not know the analytical form of the large 
deviation function, we show that violations of the Gallavotti-Cohen symmetry can be used to bound the range of 
validity of numerical results. 
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